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Abstract. Parametric excitation of a Bose-Einstein condensate (BEC) can be 
realized by periodically changing the interaction strength between the atoms. Above 
some threshold strength, this excitation modulates the condensate density. We show 
that when the condensate is trapped in a potential well of irregular shape, density 
waves can be strongly concentrated ("scarred") along the shortest periodic orbits of 
a classical particle moving within the confining potential. While single-particle wave 
functions of systems whose classical counterpart is chaotic may exhibit rich scarring 
patterns, in BEC, we show that nonlinear effects select mainly those scars that are 
locally described by stripes. Typically, these are the scars associated with self retracing 
periodic orbits that do not cross themselves in real space. Dephasing enhances this 
behavior by reducing the nonlocal effect of interference. 
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Figure 1. Steady-state density profiles of parametrically driven condensates trapped 
in a stadium billiard, as obtained from a numerical solution of Eq. ([1]) at various 
excitation frequencies |12). (a) The "horizontal" scar, (b) The "bouncing ball" scar, (c) 
The "bow-tie" scar. At the lower right side of these panels we depict the corresponding 
periodic orbit. Panel (d) represents an example of the density profile of the BEC in 
the presence of dephasing. The parameters of this figure are identical to those of panel 
(c) but with short range gaussian correlated noise added to the confining potential. 

1. Introduction 

One of the hallmarks of quantum systems whose classical counterpart exhibits chaotic 
motion is the "scar" phenomenon, i.e. the concentration of wave functions along 
the short periodic orbits of the underlying classical dynamics [Ij. Scars commonly 
appear in wave equations of quantum chaotic systems such as a single particle moving 
within a potential well whose boundary is in the form of a stadium, and otherwise 
fiat. The classical counterpart of this system is a particle specularly reflected from the 
boundary. Therefore this type of single-particle two-dimensional systems, have been 
called "billiards". Most studies of scars focus on billiards. However scars may also 
appear in many-body systems, where additional factors come into play such as the 
interparticle interactions and dephasing effects. 

In this work we analyze scars in the density excitations of Bose-Einstein condensates 
(BEC) of ultracold atoms |2], confined to a potential well of irregular shape similar to 
a billiard (see e.g. Fig. 1), and parametrically excited by a periodic change of the 
interaction strength between the atoms. Recently it has been shown, both theoretically 
[21 H] and experimentally [5], that such an excitation of BECs generates waves similar 
to surface Faraday waves in a periodically accelerated fluid. With advances in trapping 
and imaging technology, bulk parametric excitations are now feasible and the density 
variations they generate may be observed in-situ [6]. 

This system can be realized, for example, by using a rapidly moving laser beam 
or a hologram which generates a potential well [7], ^(r). This potential confines the 
condensate to a cylinder whose cross-section, in the transverse direction r = {x,y), 
may be shaped as a chaotic billiard. Additional laser beams are used to confine the 
condensate in the perpendicular direction, z, in order to form a quasi-two dimensional 
geometry [HI E]. The parametric excitation of the condensate is obtained either by 
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modulating the confining potential in the z direction [TUt ITT] , or by applying a time 
dependent magnetic field near a Feshbach resonance [T^. Both efi^ectively introduce 
a time dependent component into the interaction energy between the atoms. 

Assuming that the excitations involve a large number of atoms, the system may be 
modeled by the Gross-Pitaevskii equation, 

ih^ = (1-zr) (-^V' + V{r) -/i + ^l^rj ^ + 5g{tM'^, (1) 

where F is a dimensionless phenomenological parameter accounting for the dissipation 
|15] . m is the effective mass of the atoms, = d'^ + dy is the two-dimensional Laplacian, 
H is the chemical potential, g = denotes the bare contact interaction strength, 

where a is the s-wave scattering length and dz is the effective width of the confining 
potential in the transverse (z) direct ion [11]. Sg{t) = 2acos{2ujt) represents a small 
time dependent modulation of the contact interaction, with frequency 2uj and a small 
amplitude 2a ^ g. 

In Fig. 1 we depict a few examples of the excited condensate density, |^/'(r,t)p at 
long times where a steady state excitation is reached. These results are obtained by 
numerical solution of Eq. ([1]) in a stadium billiard potential [12]. Our numerical study 
indicates that the shortest scars, shown in panels (a) and (b), are commonly observed, 
while the "bow-tie" scar associated with a self crossing orbit, shown in panel (c), rarely 
appears, as it requires a fine tuning of the parameters. We have not been able to observe 
parametric excitation of other scars or random wave functions. This behavior, as we 
show below, results from interaction effects associated with the nonlinear term in the 
Gross-Pitaevskii equation ([T]). 

This is the plan of the paper: In section [2] we present the stability analysis of the 
problem, and show that it is described by the dissipative Mathieu's equation. Next, in 
Section [3l we focus our attention on the leading instabihty of the dissipative Mathieu 
equation, and identify the most relevant nonlinear terms of the equations of motion. We 
solve these equations in order to find a steady state solution of the problem describing 
an excitation of a single mode. In section H] we consider the case where the strength of 
excitation is sufficiently large to excite several modes, and show that interactions induce 
mode locking which promote the appearance of short scars with local density pattern in 
the form of stripes. The effect of noise and dephasing on the excitation patterns will be 
discussed in section [5l Finally, in section |6] we present a numerical study of the problem, 
and in [7] we conclude. 

2. Linear analysis 

Our analysis will be limited to the experimentally relevant regime where the Thomas- 
Fermi approximation is valid. Namely when the healing length, ^ = h/^ymfi, is much 
smaller than the size of the system, L. Under this condition, the condensate density, 
|'?/'(r, in the absence of pumping, i.e. 6g = 0, is essentially constant within the 
billiard, and equals to = fi/g. 
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When the parametric pumping is apphed weakly enough, it only modifies the phase 
of the wave function. Thus il){t) = \xjj\e'^'^^^^ where (j){t) = ^ dt'5g{t'). However, as 
the strength of the excitation increases, the uniform density solution becomes unstable, 
and a parametric resonance sets in. In order to describe the behavior of the system in 
this regime, we follow the standard Bogoliubov analysis [I6] in the presence of the time 
dependent drive and dissipation. 

Let us assume the excitation strength to be sufficiently close to the instability 
threshold, and expand the wave function in the form. 



(2) 



where Wj{t) are complex amplitudes, while Uj{r) are the eigenstates of the "unperturbed 



Hamiltonian" 



(3) 



normalized such that / d'^r\uj{r)\'^ = 1. From here on we choose these eigenstates to be 
real functions. 

Substituting Eq. 



m 



multiplying the result by Uj{r), and integrating over 
space, we obtain a set of coupled nonlinear equations for the expansion coefficients, 
Wj(t). Yet in its linearized form, this set of equations decouples, such that 

where * denotes complex conjugation, and 
M 
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The above equations reduce to the damped Mathieu equation 



fi^ + ihj cos(2wt) 



0, 



(5) 



(6) 



df^ " dt 

where rjj = Re Wj is the real part of the amplitude of the j-th mode, jj = r(/i + ej)/h 
is its effective damping rate, f2j = (ej(ej + 2/i)(l + r^))//i^, is the frequency, and 
hj = afiej /{Ti^g) is the effective pumping amplitude. The stability plane of Mathieu's 
equation is depicted in Fig. 2. The vertical axis represents the excitation strength while 
the horizontal axis is the ratio of the natural frequency of the oscillator, Qj to half of 
the pumping frequency, u. The dark tongues on this plot show the regions where the 
system is unstable, namely the values of parameters where the amplitude, Wj, grows 
exponentially in time. It is evident that the leading instability, namely the first one 
encountered as the excitation increases, is associated with the tongue where half of the 
pumping frequency is close to one the frequencies, flj. 

/^From here on we shall therefore assume that the pumping strength is sufficiently 
weak and focus on the main instability where u ^ Qj. In order to describe the behavior 




Figure 2. A schematic illustration of the stability plane of the damped Mathieu 
equation ^ at small damping coefficient F <C 1. 



of the system in this case, let us perform a Bogohubov-hke [12] pseudo-rotation which 
diagonahzes the above matrix for the case of zero pumping strength, 
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Thus, the vector {bj,b*) satisfies the equation 



d_(b, 
dt \ b* 



(7) 



(8) 



(9) 
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Assuming the pumping strength to be sufficiently close to the excitation threshold, one 
may write the amplitudes in the form bj{t) = Bj{t)e~^^^ , where Bj{t) is a slowly varying 
function with respect to the frequency u. Averaging the equations flTU]) over the rapid 
time oscillations leads to: 
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Seeking a solution where both Bj and B* behave hke e'^J*, one obtains that the 
growth exponent, Uj, is given by 

af ) = -7, ± ^h:]-{uj-uj,f. (14) 

Thus the uniform density solution becomes unstable when (Xj-^' > 0, i.e. 

h, > + (15) 

In particular, when half of the excitation frequency coincides with one of the 
eigenfrequencies of the system, u = coj, the criteria for the onset of instability is hj = jj, 
or alternatively in terms of the excitation amplitude a = hgT{^ + ej)uj/{fiej). 

3. Nonlinear analysis 

Above threshold, the system oscillates at half of the pumping frequency, and its 
amplitude grows exponentially with time. Dissipation does not saturate the growth 
of the amplitude (it only shifts the excitation threshold), and one must resort to the 
nonlinear terms of the governing equations in order to describe the evolution over long 
time scales. 

To identify the most important nonlinear terms of Eq. (|4]), notice that these may 
be divided into two categories: Quadratic terms, Wj>{t)wj"{t), proportional to integrals 
of the form / d'^ruj{r)uj'{r)uj"{T), and cubic terms, Wj'{t)wj"{t)wj"'{t), proportional 
to / d'^ruj{r)uj'{r)uj"{r)ujiii{r). At high excitation frequencies, the spatial oscillatory 
nature of the eigenstates Uj{r) implies that for any choice of the indices j, j', and j", 
the quadratic contribution is small compared to the cubic terms where either one of the 
following conditions is satisfied: (a) j = j' and j" = j'"; (b) j = j" and j' = j'"; (c) 
j = j'" and j' = j". Other cubic terms will also be small. 

Therefore, from here on, we shall keep only the leading nonlinear terms, an approach 
which is similar in spirit to the Hartree-Fock approximation. If furthermore we neglect 
the nonlinear terms that are proportional to Sg{t), and express the result in terms of 
the slowly changing amplitudes, Bj, we obtain that 

9,S, = (9i5,)iij, + (9i5,)int. (16) 

Here (c?t5j)jj^ represents the linear contribution expressed by Eqs. ( fTTl) and ( fT2l) . while 

/ dBj \ {i + r)/i ^ + 9 
k l^^n' 



dt 
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X 
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where 



Ijj' = I d'^ru^-(r)u'^,(r) 



is the nonlinear contribution which couples Bj to other modes, Bj 
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Let us consider first the steady state solution, dtBj = 0, where only one mode, say 
the j-th one, is excited. Equation f lT6|) in this case reduces to: 

[-7, + t{co - a;,)] B, - %h,e'' B] - ii±£kki|l/^.^.|5,f 5, = 0, (19) 
and its solution for = is 



^Vej(ej +2/i) 



AT. = ^^thl3±_ (20) 



/i|(l + r2) - (7, +Tuj- TiOj)^ - r7, + uj-iOj 



X 

This solution is valid above the excitation threshold (fTSjl . otherwise Nj = 0. 

Nj is roughly proportional to the number of atoms in the excited mode Uj{r). Its 
relation to the spatial variation of the BEC density follows from 

5p{v,t) = \i,{l + w,{t)u,{v)\'~\^P\' 

~ 2\ipW{r) Re Wj{t), (21) 

where 

^^'^^^ " (1 + 2/i/e-)^/^ ~ ^^^^ 

Here (pj is the phase shift between the pump and the excited wave which may be 
computed from the solution of Eq. (1191) . In particular, in the limit F — )■ it approaches 
7r/2, implying that the saturation mechanism is that of phase mismatch between the 
pump and the system. 



4. Scars 



Until now we have considered only the case where one mode is excited. However, 
usually several modes will be simultaneously excited when the system is large enough. 
The reason is that the frequency difference between neighboring modes (being inversely 
proportional to the system's area) is small compared to the width of the leading 
instability tongue near the threshold, and therefore near the minimum of the instability 
tongue (see Fig. 3) the condition (fTSll will usually be satisfied for several eigenfrequencies, 

Uj. 

For sufficiently weak excitation and a fine tuning of the frequency, it is possible to 
excite a single eigenstate of the system. However, in this case, the growth rate exponent 
is extremely small, and it takes rather long time for the system to reach a steady state. 
During this time, external noise (discussed later on) and non-parametric excitations, 
generated near the edge of the system, will overwhelm the response. 

We shall therefore consider the more practical situation where several modes are 
simultaneously excited. Our numerical simulations, as well as experimental data for 
Faraday waves in hydrodynamic systems [17], indicate that in this case the natural 
density patterns that are parametrically excited are the short scars of the system. Scars 
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Figure 3. A magnification of the leading instability tongue shown in Fig. 2, and 
the eigenfrequencies of the systems (designated by the thin black bars. When the 
excitation strength a is sufficiently high, a large number of modes are excited. The 
gray Lorenzian-like curves represent the resonances associated with some scar. 



are not exact eigenstates of the system. However they may be viewed as resonances 
whose width is determined by the typical time that the underlying classical particle 
stays in the vicinity of the periodic orbit associated with the scar [18]. 

The set of frequencies, a;f,'\ associated with a given scar (denoted by the subscript 
z/), approximately satisfies the quantization condition of a particle in a one dimensional 
square well potential whose width equals to the length of the scar, namely, wf,'-* ~ 
-^(27r/ — where c is the sound velocity, L^, is the total distance of a single round trip 
along the periodic orbit, / is an integer, and (p^ is the the Maslov phase accumulated 
along the periodic orbit (e.g. due to bouncing from the boundary of the well and passing 
through focusing points). Thus the frequency gap between two neighboring frequencies 
of short scars, Au!i, = oj^'^ — w^'"^^ = 2110/ L is, in general, much larger than the mean 
frequency spacing of the exact eigenmodes of the system, as illustrated in Fig. 3. 

In order to understand how scars may emerge from the simultaneous excitation of 
several states, it is instructive to describe the general structure of the wave functions 
of chaotic systems. These may be roughly viewed as a combination of two components: 
The first, iprand is approximately random Gaussian wave function which results from 
a superposition of plane waves in arbitrary directions and arbitrary phases[19]. This 
contribution may be associated with long orbits of the underlying classical dynamics 
and is characterized by universal statistical properties. The second component, ipscar is 
nonuniversal since it emerges from constructive interference along the shortest periodic 
orbits of the classical system. This contribution leads to the scarring patterns that 
appear in individual wave functions[lj. However, usually, scars are weak addition to 
the random background, and their effect within a single wave function diminishes as 
the semiclassical limit of short wavelength is approached. On the other hand, being a 
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resonance, the contribution of a given scar extends over a large interval of the energy as 
illustrated in Fig. 3. Thus "mode- locking" of several consequative excited states may 
result in a rather strong scarring pattern. 

To clarify the nature of mode-locking in our system, let us assume that within 
the frequency range of the excited states there is a single dominant scar, so that the 
corresponding wave functions may be written in the form [20] 

Mj(r) = aji)scar{r) + ^rand,j{r)- (23) 

Here ipscari^^) is the same scar contribution shared by all wave functions within the 
resonance and without loss of generality we shall assume that the weights of this 
contribution, a^, are real and positive. The second contribution 'iprandjij) is the 
random component mentioned above, which for different wave functions is, essentially, 
statistically independent. Consider now a superposition of these wave functions: 

u{v) = Y.q,u,{v). (24) 

i 

This superposition exhibits mode locking if the weights qj have the same phase, since in 
this situation the scar contribution adds up coherently and overwhelms the incoherent 
contribution of the statistically independent random components. 

Now, mode locking will be generated if the resulting excitation pattern is stronger 
than other possible random like patterns generated by different combinations of the 
eigenstates Mj(r). The strength of the excitation, in turn, is determined by the 
interaction matrix elements of the excited state, as follows from formula fl20|) . Namely, 
the smaller the interaction, the stronger the excitation. Thus to close our argument for 
mode locking, we have to show that the matrix elements of certain scars are smaller 
than those of random wave patterns. 

Let us consider first the interaction matrix element associated with a random wave 
function, denoted henceforth by Irand- Approximating such a wave function by a random 
superposition of plane waves moving in arbitrary directions, we obtain from formula ( |T8l) 

hand ~ 3/^, (25) 

where A is the system's area. 

Consider now the interaction matrix element associated with a stripe pattern 
(denoted by Igt) that extends, essentially, all over the system. The wave function 
associated with this pattern, u{t) ~ Y^cos(k ■ r), may be viewed as two plane waves 
moving in opposite directions. The resulting matrix element, 

/,i~3/(2^), (26) 

is half that of the random pattern. Therefore excitation of stripes, in general, is 
preferable over the excitation of random patterns. Now, both scars, the horizontal 
(Fig. la) and the bouncing ball scar (Fig. lb), being associated with self retracing 
orbits representing waves moving in opposite directions, result in a local stripe pattern. 
Therefore whenever the scar contribution, ipscar{T^)i is associated with one of these orbits, 
the system will be mode-locked in order to generate the stripe pattern. 
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Figure 4. The effective area of a scar: Tfie scar wave function is associated witli closed 
orbits that are small deformations of the central periodic orbit. Here we approximate 
such an orbit as an isosceles triangle whose base is the original horizontal periodic 
orbit. The triangle height of the longest orbit that contribute to the scar is determined 
from the requirement for constructive interference of the orbits, thus the total length 
of the longest orbit is up to one wavelength A longer than that of the periodic orbit. 
This conditions implies that that the height is V XL and therefore the effective area of 
the scar is lVXL. 



Yet, this argument is based on the assumption that the stripe pattern of the scars 
extends over the whole system. In practice the effective area of this pattern, Aeff, is 
smaller, and therefore A in Eq. f l26|) should be replaced by the proper effective area of 
the scar pattern. The bouncing ball scar occupies the central zone of the billiard whose 
area is approximately Aeff = 4:A/{4: + vr). Thus the interaction matrix element of this 
scar is 

T 3(4 + tt) 

This matrix element is still smaller than that the of the random wave function f l25|) . 
and therefore the excitation of the bouncing ball stripe pattern is stronger than that of 
a typical random pattern. Yet, the bouncing ball scar represents a special case, since 
it is associated with a continues family of periodic orbits bouncing between the two 
straight edges of the stadium billiard. In general, the effective area of a scar coming 
from a generic (isolated) periodic orbit shrinks as the excitation frequency increases. In 
particular, the effective area occupied by the horizontal scar is of order Aeff ~ V XL^ 
where A is the wavelength and L is the length of the billiard. This estimate for the 
scar's area is explained in Fig. 4. 

Replacing A in Eq. fl26|) by the effective area we obtain that the matrix element of 
the horizontal scar is 

To obtain this result we have used the relation A = (L^/16)(4 + vr) between the area 
and the length of the billiard. 

Comparing this matrix element with that of the random pattern ( l25l) it is evident 
that for sufficiently small A, i.e. high enough excitation frequency, Irand < In- Thus 
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the horizontal scar pattern is preferable when A > (4 + 7r)^/^-\/^/256 ~ L/30, or at the 
frequency range where 



where c ~ y/i/m is the sound velocity. 

What about other scars? Our numerical study indicates that scars different from 
those shown in the upper panels of Fig. 1, are rarely excited. One of these rare 
examples is the bow tie scar shown in Fig. Ic. It is somewhat difficult to associate 
this scarring pattern with the classical orbit shown in the lower right inset, however a 
Fourier transform of this wave pattern confirms that it is approximately composed of 
four components of plane waves at directions which differ by 90*^. Such a combination of 
plane waves results in square pattern and assuming it to extend over the whole billiard 
area, we obtain that the corresponding interaction matrix element is Igq = 9/(4^). It is 
larger than that associated with stripes but smaller than the matrix element of random 
wave function. Therefore, in a frequency regime which does not contain the frequencies 
of the horizontal and the bouncing-ball scars, such a scarring pattern may emerge. 
Notice, however, that similar to the horizontal scar, at high excitation frequencies the 
effective area of this scar diminishes and the random wave pattern becomes preferable. 

Longer scars are associated with rays which move in many directions and therefore, 
their interaction matrix elements become closer to that of the random wave pattern and 
thus unlikely to be observed. For instance hexagonal pattern emerging from from six 
plane waves yields an interaction matrix element Ihex = 5/{2A). 

To conclude, the horizontal and bouncing ball scars are the dominant excitation 
patterns because they generate stripes oriented perpendicular to these self retracing 
orbits. Longer periodic orbits generate more complicated density pattern due to 
interference of waves moving in several directions, e.g. at the vicinity of self crossing 
points of the periodic orbits, or near the boundary of the system where the incidence 
angle of the impinging orbit differs from zero. Near such points, the local scarring 
pattern is expected to more complicated. For example in the vicinity of a self crossing 
point where the orbit crosses itself at right angle the local density excitation will be 
in the form of a square pattern. Interactions do not favor these patterns and therefore 
scarring patterns associated with such orbits are less likely to be excited. Finally, all 
scars associated with isolated periodic orbits disappear at sufficiently high frequencies. 
The bouncing ball scar is an exception since its effective area is independent of the 
frequency. 



5. Noise and dephasing 

The above conclusions have been drawn under the assumption that the boundary of the 
billiard plays a crucial role in selecting the patterns of the excited waves in the system. 
This is because these patterns emerge from the single particle wave functions Uj{r), see 
Eq. ([3]). Yet, when the system is large enough one may expect the boundary to play a 
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minor role, and the local interactions among waves to be the dominant mechanism of 
pattern selection. 

What is the criterion of a large system? Pure dissipation in the manner included in 
Eq. ([1]) does not introduce a length scale which competes with the size of the system. 
Increasing the dissipation mainly lifts the excitation threshold but does not affect the 
role of the boundary. On the other hand, being an outcome of interference effect, scars 
are very sensitive to dephasing, namely external, random-like, forces which affect the 
phase coherence of the system. 

Dephasing emerges, for example, in the presence of noise, associated with a 
stochastic time-dependent component of the confining potential, 5l^(r,t), or interaction 
with the incoherent component of the condensate |2l]. This introduces an additional 
length scale to the problem: The dephasing length which is the typical distance 
over which a wave losses its coherence. When the dephasing length is much smaller 
than the system size, L , the BEG density patterns will be governed by local 

interactions between waves, and the role of the boundary shape will diminish as the 
system size increases. In this case, as we showed previously, the main candidates for the 
local patterns are stripes which represent the smallest interaction matrix element. 

This implies, in particular, that the horizontal and the bouncing-ball scars of the 
stadium billiard are rather stable against dephasing, while the bow tie scar shown in 
Fig. Ic is much less stable. Fig. Id shows the excitation pattern obtained by adding noise 
to the system which otherwise generate the pattern shown in Fig. Ic. Fourier transform 
confirms that this pattern is dominated by stripes. The asymmetry of this pattern is 
a clear indication that dephasing effects take place. Asymmetric excitation patterns 
appear also when adding noise to the horizontal and bouncing ball scars, however, the 
noise power needed in this case is stronger by two orders of magnitude than that required 
in order to destroy the pattern shown in Fig. Ic. 

6. Numerical study 

To check the above conclusions, we have performed an extensive numerical |24j study 
of Eq. ([1]). Some of the results of this study are shown in Fig. 1 where we considered 
a standard stadium billiard [R = 12.0 microns), with = 5 x 10^ ^"^Rb atoms, and 
assuming a strong (trapping frequency of 8.5 kHz) confinement in the z direction. The 
convergence to steady state excitation is slow, and a rather long (~ SOO/w) transient 
response of the system to the parametric drive is generic. We find an exponential 
sensitivity of the excited modes to F and a, with thresholds clearly visible. We also 
confirm that the horizontal and the bouncing ball scars patterns are the most likely 
ones to be excited. 

Dephasing has been introduced to the system by adding a time dependent 
contribution to the confining potential which is random both is space and time and 
short range correlated in both cases. We verify that as the dephasing noise increases, 
scars become less prominent, and for small size systems, the BEG density patterns do 
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Figure 5. Left panel: A billiard made of four touching circles, and the corresponding 
short periodic orbits of the system. Right panel: A typical density pattern (same 
number of atoms and density as in Fig. 1, w = 2.86/i/?i, a — 0.14^ and T = 0.03, 
reflecting the spontaneous breakdown of symmetry between the scars associated with 
the shortest periodic orbits of the system (the solid lines in the left panel). 

not respect the symmetry of the system and are dominated by stripes, as demonstrated 
in Fig. 1(d). These patterns are stable against changes in the strength of the noise, 
and its correlation length, as long as the dephasing length is smaller than the size of the 
system. 

As an additional check of the implication of our theory, we study a different 
geometry characterized by the symmetry of the dihedral point group D^^. The boundary 
of this chaotic system, formed by four circles, is shown in the left panel of Fig. 5, together 
with a few of its short periodic orbits. The symmetry of the system generates an exact 
degeneracy of some of its eigenstates. For instance, a degeneracy between pairs of states 
scarred along the two shortest periodic orbits of the system (shown by the solid lines 
in Fig. 5). This degeneracy implies that both states should be excited simultaneously. 
However, the two scars cross in the middle of the system and therefore a simultaneous 
excitation of both of them will result in a square pattern near the central point of the 
system. On the other hand, our analysis of the interaction matrix elements suggests 
that stripe patterns, associated with the excitation of a single scar, are preferable over 
the square pattern resulting form the excitation of both scars. Thus one expects to 
observe excitation patterns in which spontaneous breakdown of symmetry occurs and 
one scar is selected. This is indeed the case in many patterns parametrically excited in 
the system. An example is shown in the right panel of Fig. 5. 

Yet, the interaction matrix element is not the sole factor which determines the form 
of the excitation pattern. The precise value of the excitation frequency, for instance, also 
plays an important role, as follows from formula (l20l) . Accordingly, in some of the cases 
we observe other patterns associated with the short periodic orbits of the system, as 
shown in the upper row of Fig. 6. Our numerical simulations cannot answer the question 
whether these are metastable patterns or the the final excitation patterns. But, when 
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Figure 6. Parametric excitation patterns of BEC in the potential well shown at 
the left panel of Fig. 5. In the upper row of the figure depicted density plots 
of the excited patterns at three different excitation frequencies (from left to right, 
(jj = 3.42, 2.64, 3.7/z/?i) but the same damping and amplitude of excitation (a = 0.14g 
and F = 0.03). Below each one of these panels we show the density plot obtained for 
the same parameters but in the presence of noise which generates dephasing. 

introducing noise into the system, these patterns become unstable and taken over by 
the stripe pattern associated with one of the shortest periodic orbits, see lower row of 
Fig. 6. 




Figure 7. Parametric excitation patterns for a circular billiard. We use the same 
parameters as Fig. 6, with w = 2.7/ti/?i in (a) and (c) and to = 3.2/x/S. in (b) and (d). 
In (c) and (d) we also apply noise as in Fig. 6 

Finally, it is instiTictivc to compare these results to the case of billiard in which 
the underlying dynamics is integrable, such as a circular billiard. In Fig. 7 we show 
the results of such a simulation at two different drive frequencies (a-b), including also 
the effects of noise (c-d). We observe that the structure of the parametric excitation 
is not sensitive to noise (up until it is completely destroyed by it). This is consistent 
with our view that local stripe structure is very stable. In this billiard the stripes only 
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significantly interact at the origin. Therefore noise only causes local changes near the 
origin (see Fig. 7b and 7d). 

7. Summery 

To conclude, we show that scars in BECs result from mode-locking of many excited 
eigenstates. This effect is particularly efficient for the very short scars which exhibit 
stripe patterns. Longer scars, associated with waves moving in many directions, are 
hindered by nonlinear effects. Except for cases like the bouncing ball scar, or more 
generally scars associated with continuous families of periodic orbits, all other scars are 
expected to disappear at sufficiently high frequency (small wavelength). This implies 
that random wave patterns in EEC's may appear, in principle, only at sufficiently high 
excitation frequencies. Noise also hinders the formation of patterns other than stripes. 
The latter property may be exploited in order to measure the dephasing time in BECs. 
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